Skip to content

Conversation

Copy link

Copilot AI commented Nov 14, 2025

Pull request type

  • Code changes (bugfix, features)
  • Code maintenance (refactoring, formatting, tests)

Checklist

  • Tests for the changes have been added (if needed)
  • Docs have been reviewed and added / updated
  • Lint (black rocketpy/ tests/) has passed locally
  • All tests (pytest tests -m slow --runslow) have passed locally
  • CHANGELOG.md has been updated (if relevant)

Current behavior

Drag coefficients are limited to 1D functions of Mach number. CFD simulations typically produce drag data varying with Reynolds number and angle of attack, which cannot be properly represented.

New behavior

Drag coefficients now support N-dimensional structured grids via Function.from_grid(). Flight simulations automatically calculate and pass Reynolds number (from velocity, density, viscosity, rocket diameter) and angle of attack (from freestream velocity vectors in body frame) to multi-dimensional drag functions.

Implementation

Function Class (rocketpy/mathutils/function.py)

  • Added from_grid() class method using scipy.interpolate.RegularGridInterpolator
  • Added 'linear_grid' interpolation type with edge clamping (no extrapolation beyond bounds)
  • Grid data stored as flattened array for serialization compatibility

Flight Class (rocketpy/simulation/flight.py)

  • Added __get_drag_coefficient() helper detecting function dimensionality
  • Calculates Re = ρVD/μ and α from body-frame velocities during integration
  • Updated udot_rail1(), u_dot(), u_dot_generalized() to use helper

Rocket Class (rocketpy/rocket/rocket.py)

  • Preserves Function objects when passed as drag coefficients (previously wrapped in 1D Function)

Usage Example

import numpy as np
from rocketpy import Function, Rocket

# Create 3D drag data from CFD
mach = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
reynolds = np.array([1e5, 5e5, 1e6])
alpha = np.array([0.0, 2.0, 4.0, 6.0])  # degrees

M, Re, A = np.meshgrid(mach, reynolds, alpha, indexing='ij')
cd_data = 0.3 + 0.1*M - 1e-7*Re + 0.01*A

# Create structured grid function
drag_func = Function.from_grid(
    cd_data, 
    [mach, reynolds, alpha],
    inputs=['Mach', 'Reynolds', 'Alpha'],
    outputs='Cd'
)

# Use directly in Rocket - Flight automatically handles M, Re, α
rocket = Rocket(..., power_off_drag=drag_func, power_on_drag=drag_func)

Testing

  • 7 unit tests for Function.from_grid() covering 1D/2D/3D grids, extrapolation, validation
  • 3 integration tests for 3D drag in flight simulations
  • All 593 existing tests pass (full backward compatibility)

Breaking change

  • No

Existing 1D drag curves (Mach only) work unchanged. Multi-dimensional support is opt-in via from_grid().

Additional information

Resolves #874

CodeQL security scan: 0 alerts

Original prompt

This section details on the original issue you should resolve

<issue_title>ENH: Drag curve as a function of Re and alpha as well</issue_title>
<issue_description>Is your feature request related to a problem? Please describe.

As of now drag curve is a function of speed (or mach number). Since a lot of rocketeers use CFD simulations, it is expected that drag values won't be just a function of mach number. Drag of any body moving in fluid can depend on many major parameters. Some of them are listed below:

  • Reynolds Number (Re)
  • Angle of Attack (alpha)
  • Altitude/Air Properties
  • Rocket Geometry/Shape
  • Surface Roughness

Re accounts for the viscous effects of the flow (skin friction vs. pressure drag). Reynolds number depends on velocity, characteristic length, air viscosity, and density, and is crucial for modeling drag at low speeds and during laminar-turbulent transitions. Non-zero angles of attack change the effective cross-sectional area and cause flow separation or vortex formation, increasing drag. Altitude or atmospheric variation can be handled by RocketPy separately to update Re or atmospheric data rather than being built directly into the drag curve.
Drag could be input as Cd(M, α) or Cd(M, Re, α).

Describe the solution you'd like

  1. A multidimensional array implementation of Cd=f(M,Re,α) in Rocket Class.
  2. Extension of drag_coeff in flight class to allow multiple params as input - M, Re and alpha. Probably implement new get_value_opt like function specifically for this.
  3. Check if current interpolation and extrapolation methods of the Function class are applicable in case of 2D/3D arrays. If not then update the same.
  4. Ensure backwards compatibility to still allow for use of Cd=f(M).

Additional context

It does not make sense that the drag curve will be given as an input of Cd=f(Re) or Cd=f(α). Physically these are not capturing the entire flow physics at the same accuracy at which Cd=f(M) does as mach number based effects like shock wave propagation become more important in rocket flights. May be still allow for single array functionality with Re and α but provide a note that the results may not be accurate!</issue_description>

<agent_instructions>Here is a specific, technical prompt you can use to delegate this issue to the GitHub Copilot agent. This prompt directs it to the correct implementation strategy we discussed (enhancing the Function class).


Optional Prompt for GitHub Copilot

Here is a specific prompt to guide the implementation for Issue #874.

1. Core Objective:
Implement the enhancement described in issue #874. The goal is to allow the drag coefficient Cd to be defined as a 3D function of Mach (M), Reynolds number (Re), and Angle of Attack ($\alpha$).

2. Implementation Location:
The primary changes must be made within the rocketpy.Function class (in rocketpy/function.py). This is not just a change to the Flight or Rocket class; the Function class itself must be enhanced to handle N-dimensional structured grid data.

3. Add New Interpolation Method:
The current shepard interpolation is inefficient for structured CFD grids.

  • Add a new interpolation option (e.g., interpolation='linear_grid' or interpolation='grid') to the Function class __init__.
  • This new method MUST use scipy.interpolate.RegularGridInterpolator as its backend. You will need to import this from scipy.interpolate.
  • The Function class state should store the initialized RegularGridInterpolator object to be called later.

4. Add New Function.from_grid() Class Method:
Manually creating a flattened list of points is impractical for users.

  • Create a new class method: @classmethod def from_grid(cls, grid_data, axes, inputs, output, interpolation='linear_grid', **kwargs).
  • This method will accept:
    • grid_data: An N-dimensional NumPy array (e.g., a 3D array for $C_d$).
    • axes: A list of 1D arrays representing the axes of the grid (e.g., [M_axis, Re_axis, A_axis]).
    • inputs: A list of strings, e.g., ['Mach', 'Reynolds', 'Alpha'].
    • output: A string, e.g., 'Cd'.
  • Inside this method, you will initialize the scipy.interpolate.RegularGridInterpolator with the provided axes and grid_data and store it in the new Function object.

5. Handle Extrapolation (Crucial):
The new grid interpolation must clamp to the edge (i.e., use the nearest edge value) for out-of-bounds requests. It should not extrapolate.

  • When initializing RegularGridInterpolator, set bounds_error=False.
  • The fill_value should be set to None, which (when bounds_error=False) causes the interpolator to use the edge values for out-of-bounds points.

6. Update Flight Class:
The Flight class (likely in flight.py or related files) must be modified to use this new function.

  • The flight simulation loop needs to detect the type of drag function it has.
  • If the `rocket.aerodynamic_surfaces.dra...

💬 We'd love your input! Share your thoughts on Copilot coding agent in our 2 minute survey.

@Gui-FernandesBR Gui-FernandesBR changed the base branch from master to develop November 14, 2025 19:33
Copilot AI changed the title [WIP] Enhance drag curve to include Re and alpha variables ENH: Add multi-dimensional drag coefficient support (Cd as function of M, Re, α) Nov 14, 2025
Copilot finished work on behalf of Gui-FernandesBR November 14, 2025 20:04
@aZira371
Copy link
Collaborator

aZira371 commented Nov 15, 2025

The code implementation is as I planned it to be done. There might be a slight edge case which should be addressed in regular grid interpolator unit tests or with warning: that the grid was not sorted in ascending order as intended. This can happen with repeated values in the grid. Apart from that the changes are okay to be merged! @Gui-FernandesBR @copilot

@Gui-FernandesBR Gui-FernandesBR force-pushed the copilot/enhance-drag-curve-functionality branch from 356e14a to b621a16 Compare November 15, 2025 16:20
@aZira371 aZira371 self-assigned this Nov 15, 2025
@aZira371
Copy link
Collaborator

@Gui-FernandesBR Is there a way to see this branch on VSC? I am unable to access it there. It is not being listed as a remote branch. I think one of the previous lint errors was related to the new grid method being too big in function.py. It should be ideally broken down into helper functions or pylint has to be disabled for it.

@codecov
Copy link

codecov bot commented Nov 24, 2025

Codecov Report

❌ Patch coverage is 89.03226% with 17 lines in your changes missing coverage. Please review.
✅ Project coverage is 80.36%. Comparing base (9cf3dd4) to head (74fe825).
⚠️ Report is 6 commits behind head on develop.

Files with missing lines Patch % Lines
rocketpy/mathutils/function.py 82.79% 16 Missing ⚠️
rocketpy/simulation/flight.py 98.21% 1 Missing ⚠️
Additional details and impacted files
@@             Coverage Diff             @@
##           develop     #875      +/-   ##
===========================================
+ Coverage    80.27%   80.36%   +0.08%     
===========================================
  Files          104      106       +2     
  Lines        12769    13122     +353     
===========================================
+ Hits         10250    10545     +295     
- Misses        2519     2577      +58     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@aZira371
Copy link
Collaborator

The pull request is ready for review now. Pylint is failing due to a problem in .pylintrc, I am not sure if this file should be corrected in this PR or not as this is an unrelated error. Tests did pass locally on my end before merging from develop. @RocketPy-Team/code-owners will need help from you, because different tests are failing at different places - on github test_sensor failed but on local after merging develop test_environment failed.

@aZira371 aZira371 marked this pull request as ready for review November 24, 2025 20:31
@aZira371 aZira371 requested a review from a team as a code owner November 24, 2025 20:31
@aZira371 aZira371 requested a review from phmbressan November 27, 2025 05:09
Copilot AI and others added 11 commits November 27, 2025 09:24
Co-authored-by: Gui-FernandesBR <[email protected]>
- MNT: velocity_body was not being used in get_drag_coefficient, removed it as an input
… test_multidim_drag.py

-MNT:  removed unused velocity in body frame parameter requirement from all instances of get_drag_coefficient in flight.py

- MNT: corrected docstring for get_value_opt_grid in function.py

- MNT: shifted import of classes before the definition of functions in test_multidim_drag.py
- MNT: rearranged the docstring of from_grid in function.py to match the expected output of doctest
- MNT: reran make format and lint on function.py to correct after previous changes to from_grid
- MNT: disables pylint unused private member for get_value_opt_grid as it is called upon dynamically by from_grid

- MNT: disabled pylint too many statement for from_grid for now and added a to-do to refactor it into smaller methods/helper functions

- MNT: updated .pylintrc to record Re as good name
- MNT: Re variable was unused in test_3d_drag_with_varying_alpha thus replaced it
* Add tests for shepard_fallback in test_function_grid.py

Co-authored-by: aZira371 <[email protected]>

---------

Co-authored-by: copilot-swe-agent[bot] <[email protected]>
Co-authored-by: aZira371 <[email protected]>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@aZira371 is linear_grid really the best name here? regular_grid should be more similar to scipy, no?

Copilot finished reviewing on behalf of Gui-FernandesBR November 27, 2025 12:30
Copy link

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR implements multi-dimensional drag coefficient support, allowing drag to vary with Mach number, Reynolds number, and angle of attack. This enhancement enables RocketPy users to utilize CFD simulation data that naturally depends on multiple flow parameters, significantly improving aerodynamic modeling accuracy.

Key changes:

  • Added Function.from_grid() class method for creating N-dimensional structured grid interpolators using scipy.interpolate.RegularGridInterpolator
  • Modified Flight class to automatically calculate Reynolds number and angle of attack during trajectory integration, passing these to multi-dimensional drag functions
  • Updated Rocket class to preserve Function objects when provided directly as drag coefficients (backward compatible with existing 1D usage)

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 7 comments.

Show a summary per file
File Description
rocketpy/mathutils/function.py Adds from_grid() class method, linear_grid interpolation type, and __get_value_opt_grid() for structured grid evaluation with edge-clamping extrapolation
rocketpy/simulation/flight.py Implements __get_drag_coefficient() helper to detect function dimensionality and calculate Re/α during integration; updates all u_dot methods
rocketpy/rocket/rocket.py Preserves Function objects when passed as drag coefficients instead of wrapping them, enabling multi-dimensional drag
tests/unit/mathutils/test_function_grid.py Comprehensive unit tests for grid interpolation (1D/2D/3D), extrapolation modes, validation, and fallback behavior
tests/integration/test_multidim_drag.py Integration tests verifying 3D drag works in full flight simulations and maintains backward compatibility with 1D drag
.pylintrc Adds 'Re' (Reynolds number) to accepted variable names

Comment on lines +144 to +149
def test_shepard_fallback_warning():
"""Test that shepard_fallback is triggered and emits a warning.
When linear_grid interpolation is set but no grid interpolator is available,
the Function class should fall back to shepard interpolation and emit a warning.
"""
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The docstring for test_shepard_fallback_warning states that "the Function class should fall back to shepard interpolation and emit a warning," but the test implementation doesn't actually test the fallback scenario correctly. The test manually changes interpolation to 'linear_grid' on a Function that wasn't created with from_grid(), which is not a realistic scenario.

Consider renaming the test to better reflect what it actually tests (e.g., test_linear_grid_without_grid_interpolator_warns) or restructure it to test a more realistic fallback scenario if one exists.

Copilot uses AI. Check for mistakes.
mu = self.env.dynamic_viscosity.get_value_opt(z)
freestream_speed = np.linalg.norm(freestream_velocity_body)
characteristic_length = 2 * self.rocket.radius # Diameter
reynolds = rho * freestream_speed * characteristic_length / mu
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Potential division by zero: reynolds = rho * freestream_speed * characteristic_length / mu on line 1412 doesn't check if mu (dynamic viscosity) is zero or near-zero before division. While atmospheric models typically provide non-zero viscosity values, defensive programming suggests adding a check or using a safe division operation to prevent crashes in edge cases.

Consider: reynolds = rho * freestream_speed * characteristic_length / max(mu, 1e-10) or adding an explicit check.

Suggested change
reynolds = rho * freestream_speed * characteristic_length / mu
# Defensive: avoid division by zero in case mu is zero or near-zero
reynolds = rho * freestream_speed * characteristic_length / max(mu, 1e-10)

Copilot uses AI. Check for mistakes.
Comment on lines +100 to +126
def test_3d_drag_with_varying_alpha():
"""Test that 3D drag responds to angle of attack changes."""
# Create drag function with strong alpha dependency
mach = np.array([0.0, 0.5, 1.0, 1.5])
reynolds = np.array([1e5, 1e6])
alpha = np.array([0.0, 5.0, 10.0, 15.0])

M, _, A = np.meshgrid(mach, reynolds, alpha, indexing="ij")
# Strong alpha dependency: Cd increases significantly with alpha
cd_data = 0.3 + 0.05 * M + 0.03 * A
cd_data = np.clip(cd_data, 0.2, 2.0)

drag_func = Function.from_grid(
cd_data,
[mach, reynolds, alpha],
inputs=["Mach", "Reynolds", "Alpha"],
outputs="Cd",
)

# Test at different angles of attack
# At zero alpha, Cd should be lower
cd_0 = drag_func(0.8, 5e5, 0.0)
cd_10 = drag_func(0.8, 5e5, 10.0)

# Cd should increase with alpha
assert cd_10 > cd_0
assert cd_10 - cd_0 > 0.2 # Should show significant difference
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test test_3d_drag_with_varying_alpha verifies that the drag function responds to different alpha values when called directly, but doesn't actually verify that the flight simulation is using angle of attack correctly during integration. Consider adding assertions that check the drag coefficient values during an actual flight simulation (e.g., by examining flight data or verifying that a rocket with high-alpha-sensitivity drag behaves differently than one without).

Copilot uses AI. Check for mistakes.
Comment on lines +4085 to +4088
extrapolation : str, optional
Extrapolation behavior. Default is 'constant', which clamps to edge values.
'constant': Use nearest edge value for out-of-bounds points.
'zero': Return zero for out-of-bounds points.
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The docstring states that extrapolation can be 'constant' or 'zero', but the __get_value_opt_grid method (line 726-727) allows other extrapolation types (like 'natural') to fall through to the interpolator's default behavior. The docstring should document what happens with 'natural' extrapolation (linear extrapolation from edge gradients) or explicitly state that only 'constant' and 'zero' are supported for grid interpolation.

Additionally, consider adding validation in from_grid() to warn or error if an unsupported extrapolation type is specified.

Copilot uses AI. Check for mistakes.
Comment on lines +4207 to +4209
func.y_array = func._image[: len(axes[0])] # Placeholder
func.y_initial, func.y_final = func._image[0], func._image[-1]
if len(axes) > 2:
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Potential issue with y_array initialization for multi-dimensional grids. Line 4207 sets func.y_array = func._image[: len(axes[0])] which only takes the first n values from the flattened grid where n is the length of the first axis. This creates an inconsistent y_array that doesn't properly represent the function data.

For grid functions with dimensions > 1, y_array is not a meaningful concept (it's designed for 1D functions). Consider either:

  1. Setting y_array to the full flattened image: func.y_array = func._image
  2. Setting it to None or an empty array with a comment explaining that y_array is not applicable for grid functions
  3. Documenting this limitation clearly

This could cause issues if other parts of the code expect y_array to contain all function values.

Suggested change
func.y_array = func._image[: len(axes[0])] # Placeholder
func.y_initial, func.y_final = func._image[0], func._image[-1]
if len(axes) > 2:
# For 1D functions, y_array contains the function values.

Copilot uses AI. Check for mistakes.
Comment on lines +4185 to +4191
mesh = np.meshgrid(*axes, indexing="ij")
domain_points = np.column_stack([m.ravel() for m in mesh])
func._domain = domain_points
func._image = grid_data.ravel()

# Set source as flattened data array (for compatibility with serialization, etc.)
func.source = np.column_stack([domain_points, func._image])
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Performance consideration: Creating the full flattened domain and image arrays (lines 4185-4191) for compatibility purposes requires O(n^d) memory where n is the typical axis length and d is the number of dimensions. For large 3D grids (e.g., 100x100x100 points), this creates redundant data since the grid interpolator already has the structured data.

While this is necessary for backward compatibility with code expecting _domain and _image attributes, consider documenting this memory overhead and potentially adding a parameter to skip this flattening if compatibility with existing Function methods is not required.

Copilot uses AI. Check for mistakes.
Comment on lines +1432 to +1445
# Determine which parameters to pass based on input names
input_names = [name.lower() for name in drag_function.__inputs__]
args = []

for name in input_names:
if "mach" in name or name == "m":
args.append(mach)
elif "reynolds" in name or name == "re":
args.append(reynolds)
elif "alpha" in name or name == "a" or "attack" in name:
args.append(alpha_deg)
else:
# Unknown parameter, default to mach
args.append(mach)
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The input parameter name matching uses lowercase comparison and substring matching (e.g., "mach" in name, "reynolds" in name). While functional, this could match unintended parameter names (e.g., "machining" would match "mach").

Consider using more precise matching such as:

  • Exact matches: name in ['mach', 'm']
  • Or document the substring matching behavior clearly in the method docstring

This is a minor robustness improvement for edge cases.

Copilot uses AI. Check for mistakes.
Drag coefficient value
"""
# Check if drag function is multi-dimensional by examining its inputs
if hasattr(drag_function, "__inputs__") and len(drag_function.__inputs__) > 1:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not confortable with this comparision here.

        if hasattr(drag_function, "__inputs__") and len(drag_function.__inputs__) > 1:

This could be done within the Function class.

An attribute called is_multidimensioanl: bool should be really useful.

Then the comparision here would be:

if isinstance(foo, Function) and foo.is_multidimensional:

muche easier to read!

Copy link
Member

@Gui-FernandesBR Gui-FernandesBR left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • @MateusStano please review changes in the flight.py file
  • @phmbressan please review changes in the function.py ile

Everything is looking good, @aZira371 great work here! I'm confortable with what has been done so far. I believe we can include this one on the next release.


@aZira371 could you please update changelog?

@aZira371
Copy link
Collaborator

Thanks for the review @Gui-FernandesBR will soon work on this and update the changelog.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

Status: Backlog

Development

Successfully merging this pull request may close these issues.

ENH: Drag curve as a function of Re and alpha as well

3 participants